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We study generalizations of the singlet-sector amplitude-product (AP) states in the valence-bond 
basis of S = 1/2 quantum spin systems. In the standard AP states, the weight of a tiling of 
the system into valence bonds (singlets of two spins) is a product of amplitudes depending on 
the length of the bonds. We here introduce correlated AP (CAP) states, in which the amplitude 
product is further multiplied by factors depending on two bonds connected to a pair of sites (here 
nearest neighbors). While the standard AP states can describe a phase transition between an 
antiferromagnetic (Neel) state and a valence-bond solid (VBS) in one dimension (which we also 
study here), in two dimensions it cannot describe VBS order. With the CAP states, Neel- VBS 
transitions are realized as a function of some parameter describing the bond correlations. We here 
study such phase transitions of CAP wave-functions on the square lattice. We find examples of 
direct first-order Neel- VBS transitions, as well as cases where there is an extended U(l) spin liquid 
phase intervening between the Neel and VBS states. In the latter case the transitions are continuous 
and we extract critical exponents and address the issue of a possible emergent U(l) symmetry in the 
near-critical VBS. We also consider variationally optimized CAP states for the standard Heisenberg 
model in one and two dimensions and the J-Q model in two dimensions, with the latter including 
four-spin interactions (Q) in addition to the Heisenberg exchange (J) and harboring VBS order 
for large Q/ J. The optimized CAP states lead to significantly lower variational energies than the 
simple AP states for these models. 

PACS numbers: 75.10.Kt, 75.10.Jm, 75.40.Mg, 75.40.Cx 



I. INTRODUCTION 

The valcncc-bond (VB) basisi - — is ideally suited for 
describing many different types of ground states and 
low-energy excitations of quantum spin modelsi 7 - - — In 
the case of S = 1/2 spins in the singlet sector, a ba- 
sis state corresponds to a tiling of the lattice into bonds 
connecting pairs of sites forming singlets, such that each 
spin belongs to one bond. This basis is overcomplcte if 
bonds of all lengths are included. To describe the ground 
state of a Hamiltonian with bipartite interactions, only 
bonds connecting sites on different sublattices have to 
be included — this restricted VB basis exactly reproduces 
Marshall's sign rulei£ for the ground state of such a sys- 
tem. Thus, in this basis the wave function is positive 
definite and can be sampled using Monte Carlo (MC) 
techniques. We here investigate a class of bipartite corre- 
lated VB wave functions which can exhibit valcncc-bond- 
solid (VBS) order and related interesting quantum phase 
transitions in one and two dimensions. 

In this introductory section we provide some further 
background and motivation for studying VB states. We 
review the definition and properties of the well-studied 
Liang-Doucot-Andcrson amplitude-product states^ and 
introduce their more versatile generalizations — the cor- 
related AP (CAP) states that we focus on in this paper. 
We discuss reasons to study such states in the context 
of quantum phase transitions from the antiferromagneti- 
cally ordered Neel state into non-magnetic VBS and spin 
liquid states. 



A. Valence-bond states and Marshall's sign rule 

While some analytical work has been carried out in the 
VB basis^^ in most quantitative calculations MC sam- 
pling of the bonds must normally be used to reliably eval- 
uate expectation values. Since the basis is overcomplcte, 
the non-negative definiteness of the wave function is a 
requirement to avoid problems due to negative sampling 
weights (the sign problem). Thus, in most cases VB MC 
calculations are restricted to bipartite (non-frustrated) 
systems. A two-spin singlet (VB) connecting sites a and 
b on sublattices A and B is then defined according to the 
following phase convention: 

(a, 6) = (talfc -Iatfc)/V2. (1) 

Marshall's sign rule is then incorporated for any tiling of 
an even number N of spins into N/2 singlets, 

\V) = \(a 1 ,b 1 )---(a N/2 ,b N/2 )}, (2) 

i.e., when expressed in the standard basis of z spin com- 
poncnts, \Z) = \Sf,...,S z N ), 

l y H^X>v-(^>, (3) 

the sign of a non-zero coefficient ipv{Z), i.e., for states 
with antiparallel spins on each bond, is given by 

sigi# v (Z)] = MZ) = (-1)™ A S (4) 



where riAi is the number of 4- spins on sublattice A. The 
wave function tpo(V) of the ground state of such a system 
expressed in the VB basis, 

i*o>=5>ooom (5) 

V 

is therefore non-negative. When using MC simulations, 
e.g., with a variational wave function \^>) approximating 
|*o), this is essential, because the ovcrcomplctcncss im- 
plies that the sampling weigh is not |^(1^)| 2 , as it would 
be in a standard orthogonal basis, but ip(V)i/j(V')(V'\V) 
for simultaneously sampled non-orthogonal bra and ket 
configurations \ V) and (V'\. The overlap (V'|V) and ma- 
trix elements of operators of interest in characterizing the 
states can be easily calculated,— ~— as we will discuss later 
below. 



B. Amplitude-product states 

The most commonly used variational states in this con- 
text are the AP states introduced by Liang, Doucot, 
and Anderson.— Here one associates a bond connecting 
two sites (a, b) with an amplitude h(a, b), which in the 
case of a translationally-invariant system is only a func- 
tion of the lattice vector r a b separating the two sites; 
h(a, b) = h{r a b). The wave function coefficient for a VB 
configuration V is then 

<P(V) = l[h(rr*, (6) 

r 

where n r is the number of bonds of "shape" r in the 
configuration. 

The amplitudes h(r) can be used as variational pa- 
rameters. In the original work with the AP states to 
describe the ground state of the two-dimensional (2D) 
Heisenberg model^ only the amplitudes for a small num- 
ber of short bonds were optimized, and different func- 
tional forms (exponentially or power-law decaying with 
the distance r) were tested. In later work all the am- 
plitudes (on finite lattices) were optimized, leading to 
relative energy errors (deviations from results from unbi- 
ased quantum Monte Carlo, QMC, calculations) of less 
than 0.1%i 19 i 20 In the optimal state the amplitudes de- 
cay asymptotically as 1/r 3 , which is also the result of a 
mean-field VB approach^ 

In some cases, if one is just interested in the proper- 
ties of some class of states without reference to a specific 
Hamiltonian, the optimization step is not needed. This 
approach has been taken in recent studies of the proto- 
typical resonating VB (RVB) spin-liquid state consisting 
of the superposition of all configurations of the shortest 
(nearest-neighbor) bonds on the square lattice^ - — and 
also in the presence of some fraction of the second bipar- 
tite bond (fourth- neighbor) .— These wave functions, for 
which the parent hamiltonian was recently identified (in 
the case of nearest- neighbor bonds only) j24 has exponen- 
tially decaying spin correlations but power-law decaying 
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FIG. 1: (Color online) Configurations of short valence bonds 
(shown in red) connected to a lattice link b (indicated by 
thick black bars). Their associated CAP weights, Eq. (0, are 
Cb(ri, V2) with n = T2 = 1 and b here being a horizontal link. 
We will later use a notation with weights d for CAP states 
including only short-bond correlations, with i corresponding 
to (ri,r2) according to the labeling above. 



VBS correlations. A phase transition from the Neel state 
into this kind of spin liquid can be achieved by using am- 
plitudes of the form h(r) cx l/r K and tuning the exponent 
k to a critical value d 2 ' 25 



C. AP states with bond correlations (CAPs) 

One of the motivations of the work reported in the 
present paper is to obtain a variational description of the 
2D Neel- VBS transition. For this, we need a class of wave 
functions beyond the AP states, since they do not exhibit 
VBS order (while the ID variants do, as we will discuss in 
Sec. IIIip . The 2D non-magnetic AP states are beleived to 
always be spin liquids, with exponentially decaying spin 
correlations and power-law VBS correlations, similar the 
prototypical short-bond RVB states. 12 ' 21 

We study a class of generalized AP states defined 
with bond-correlation factors multiplying the AP wave 
function ©. We take these factors to be of the form 
(7ft[ri(6), T2(b)), where b denotes a nearest-neighbor link 
on a ID chain or 2D square lattice (or, more generally, 
any lattice with some imposed bipartition), and Ti(b), 
r 2 (6) are the shapes of the two VBs connected to this 
bond (with the case of there being just a single bond 
connecting the two sites being a special case). Thus, the 
wave-function coefficient is 

ij(V) = l[h(rr*l[C b [r 1 (b),r 2 (b)}. (?) 

r b 

For a translationally invariant system C(,(ri, r 2 ) for given 
(ri,r 2 ) depends only on the orientation (horizontal or 
vertical in two dimensions) of the bond b, and these 
weights also should obey applicable lattice symmetries. 
The number of different correlation factors is then cx iV 2 
for a system of N spins. For simplicity of the notation 
we hereafter suppress the subscript b. 

In principle, in variational CAP calculations all cor- 
relation factors can be optimized, along with the ampli- 
tudes h(r), but one can also opt to consider only those 
factors C(ri, r 2 ) for which r\, r 2 < r max , with some max- 
imum bond length r max , and set the remaining weights 
to unity. The possible two-bond configurations with 
r mM — 1 arc illustrated in Fig. [TJ In variational cal- 
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dilations one would expect the energy to decrease mono- 
tonically with r max , which we will demonstrate explicitly 
in Sec. EJ 

Beyond improving the energy in variational calcula- 
tions, the correlation factors also play an important qual- 
itative role in 2D systems — without bond correlations, 
the standard AP states are either long-range Neel or- 
dered (although they do not, by construction, break the 
spin-rotation symmetry, they can still develop the magni- 
tude of the sublattice magnetization) or arc RVB spin liq- 
uids with critical VBS correlations (as discussed above in 
the context of the short-bond RVBs). They cannot form 
VBS order. In contrast, the trivial ID AP state with 
only short bonds is an extreme case of a two-fold degen- 
erate VBS state with alternating links with or without a 
VB. This kind of long-range order remains stable also in 
the presence of some fraction of longer bonds, as we will 
discuss below in Sec. IIIII The generalized CAP states (J7J 
can exhibit 2D VBS order if the correlation factors favor 
such correlations strongly enough. This is true even with 
correlations only involving only the shortest bonds on 
the square lattice, illustrated in Fig. [TJ The CAP states 
open the possibility to study the Neel-VBS transition 
in classes of wave functions with the bond correlations 
parametrized in some way, and to carry out improved 
variational calculations for systems with VBS order or 
Neel states with significant VBS fluctuations. 



D. Purpose and outline of the paper 

One of our main reasons to study the CAP states 
here is to investigate their abilities to describe the 2D 
Neel-VBS transition. This transition has received con- 
siderable interest recently in the context of "decon- 
fined" quantum-criticality (DQC) . 26 ' 27 Following earlier 
work on VBS states and quantum-critical phenomena in 
antifcrromagnet o 10 ' 28 " — and topological aspects of phase 
transitions in 3D classical spin systems^ Senthil et 
al. propose d 26 ! 27 that the 2D Neel-VBS transition is of 
an unusual kind where the standard Landau rule stip- 
ulating a gencrically first-order transition between the 
two ordered states is violated. The two order parame- 
ters in the DQC scenario are both consequences of the 
same underlying more fundamental objects — spinons in- 
teracting with an emergent U(l) symmetric gauge field. 
The spinons condense in the Neel state and confine as 
pairs in the VBS state. Unbiased QMC studies of J- 
Q models,— ~— which include certain multi-spin inter- 
actions (Q) in addition to the standard Hcisenberg ex- 
change (J), are in general in good agreement with the 
theoretical predictions. Among the most interesting fea- 
tures observed is an emergent U(l) symmetry of the 
VBS order parameter (presumably reflecting the emer- 
gent gauge field — the "photon") as the critical point is 
approached from the VBS side. Moreover, studies of 
SU(iV) generalizations of the J-Q model^ 6 - and other spin 
models^ have allowed direct connections with analytical 



large-iV calculations for the non-compact CP N 1 field 
theory argue d 26 ' 27 to describe the transition. 



1. Scope of the paper 

We here investigate whether the 2D Neel-VBS tran- 
sition can be correctly captured with a simple ansatz 
wave-function of the form ([7} with fixed single-bond am- 
plitudes (of a power-law form) and continuously varying 
short-bond correlation weights of the form in Fig. [TJ The 
result so far is negative, in the sense that we do not ob- 
serve the same kind of continuous VBS transition as in 
the J-Q model. Instead, with parameters chosen such 
that there is direct Neel-VBS transition, we find strong 
discontinuities. In other cases we find an RVB spin liq- 
uid intervening between the ordered phases. Thus, it 
still remains an interesting challenge to find a simple VB 
description of the DQC point. 

Looking at the VBS order-parameter distribution, we 
do not observe any emergent U(l) symmetry of the VBS 
at the continuous VBS to RVB transition, impliying that 
this is not a "deconfined" transition in this case. Never- 
theless, we find interesting scaling properties of the an- 
gular VBS fluctuations, although the length-scale associ- 
ated with them are not divergent. 

In addition to the 2D studies, we also closely examine 
the Neel-VBS transition within the standard AP states 
in one dimension, with amplitudes h(r) of the form l/r a . 
This transition, which occurs at a critical value of a 
(which is not universal but depends on the detailed form 
of the amplitudes for small r) was previously studied by 
Bcach^ but only the spin correlations were computed. 
Here we extract also the VBS correlations and confirm 
that there is a single critical point versus a. The expo- 
nents are continuously varying, depending on the short- 
bond amplitudes. 

We also report variational calculations including opti- 
mization with the CAP states, minimizing the energy for 
ID and 2D Heisenberg and J-Q models. Naturally, bond 
correlations have a significant improving effect in VBS 
phases, but they help also to improve the 2D Neel state 
and the critical ground state of the ID Heisenberg chain. 



2. Outline of the paper 

In Sec. |H] we briefly describe the technical aspects of 
MC calculations with AP and CAP states. In Sec. IIIII we 
study the Neel-VBS transition in ID AP states, and in 
Sec. II VI we study the more rich set of states and quantum 
phase transitions in 2D CAP states. In Sec.|V]we present 
some ID and 2D AP and CAP variational calculations 
(minimizing the energy as a function of the amplitudes 
and correlation factors) for prototypical model Hamilto- 
nians with Neel and VBS order (the Heisenberg and J- 
Q models), showing how bond correlations improve the 
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FIG. 2: (Color online) Transition graph on a 4 x 4 lattice 
consisting of two states, V a and Vg, which are depicted by 
red and black bonds, respectively. This transition graph has 
two loops formed by alternating bonds of V a and Vp . 



states. We conclude in Sec. I VII with a summary and dis- 
cussion of future prospects. 



II. MC SAMPLING OF CAP STATES AND 
CALCULATION OF OBSERVABLES 

In an AP or CAP state with the wave function of the 
form Q the expectation value of an observable O can be 
written for the purpose of importance sampling as 



mom 



E a p w afS o a p 



(8) 



E a p w a p 

where the configuration weight is 

W ct p = ^(Vp)^(V ol )(Vp\V a ), (9) 
and the estimator corresponding to O given by 

mo\v a ) 



a, 



(Vp\V a ) 



(10) 



Here the overlap (yg|V Q ) is evaluated by counting the 
number l a p of loops in the transition graph of V a and 
Vp^ 



{Vp\V a )=2 l °">. 



(11) 



An example with two loops is shown in Fig. [5J Below 
we briefly discuss MC sampling of the VB configurations 
and estimators for some important observables. 



A. Bond and spin sampling schemes 

For a given functional form for h(r) and bond corre- 
lation factors C(ri,r2), an expectation value (O) can be 
computed stochastically by importance sampling accord- 
ing to the weight W a p- Using some random reconfigu- 
ration of bonds in either the state V a or Vp or both of 



them, the standard Metropolis acceptance probability for 
a modified configuration (V a >, Vp>) is 



-^accept — I 

where the weight ratio is 



w a , p . 

W a p 



.1 



TP, 



a'P' 



W a p 



i>(V a )i>(Vp) 



(12) 



(13) 



Even for the simplest update involving changes in just 
two bonds,— the calculation of the change in the new num- 
ber of loops (l a ' pi —l a p) can require a computational time 
up to oc N for each new update proposal, since a Neel 
state has extensive loops (while magnetically disordered 
states have only short loops). Since the number of such 
updates in each MC sweeps should also be proportional 
to N, this type of update leads to a total computational 
time 0(N 2 ) for a full sweep in a Neel state, while in a 
non-magnetic state the scaling is O(N). 

The unfavorable scaling in the Neel state can be 
avoided by working in a combined space of both spins 
and bonds^ where the VBs are also sampled, by ran- 
domly selecting cither y a 4-6 or lotb f° r eacn singlet (a, b). 
Since the spin basis is orthogonal, all spins in the bra and 
the kct have to be the same, and a consistent assignment 
for both V a and Vp, thus, implies that the spins on each 
loop in the transition graph follow a staggered, tltl 
pattern. The overlap (fTTj) in the pure VB basis then 
follows, since there are two possible staggered configura- 
tions on each loop. The spins are periodically updated 
by flipping all spins in randomly selected loops. By such 
spin sampling, the weighting by the number of loops is 
accounted for automatically, due to the entropic effect 
of favoring configurations with large numbers of loops, 
without the need for actually counting the loops. For a 
detailed description of the combined _spin-bond basis and 
simulations in it, we refer to Refs. [lj| and [2l|. Here we 
just note two different ways of updating the bond config- 
urations: 

(i) In the two-bond reconfiguration scheme, an elemen- 
tary MC move consists of choosing two sites on the same 
sublattice at random (typically the two bonds on a ran- 
domly chosen pair of next- nearest- neighbor spins), and 
exchanging the bonds connected to these two sites for 
the other possible bipartite configuration. Such a recon- 
figuration is only possible if the spin states on the two 
selected sites are the same. If that is the case, the accep- 
tance probability (|12[) is applied, where in the combined 
spin-bond basis the weight is 



W a p = TP(Vp)ll>(Va), 



(14) 



instead of Eq. (j9]), and with W a 'p' /W a p evaluated using 
only the bond amplitudes and correlation factors in ([7]) 
affected by the change. 

(ii) In a loop update, we start by removing a dimcr 
randomly from two connected sites, creating two defects 
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("holes"). We keep one defect stationary and move the 
second one by connecting one end (the one on the same 
sublattice) of a chosen bond to it (hence moving the hole 
to the previous location of that end of the bond). The 
bond to move should be chosen probabilistically in such 
a way as to satisfy detailed balance, which is relatively 
straight-forward in the case of AP state s 19 i 40 but more 
complicated when bond correlations are included. In the 
present work we have used loop updates only for pure AP 
states, while we use two-bond updates for CAP states. 
The latter are also efficient enough to study relatively 
large lattices (with thousands of spins). 

It should be noted here that VB configurations can be 
classified according to topological "winding numbers" M 
In AP or CAP states defined with only short bonds, 
the two-bond update conserves the winding number, but 
with no restriction on the bond-length such updates can 
change the winding number. In practice, if the bond 
probability (which depends on the single-bond ampli- 
tudes as well as the correlation factors in CAP states) 
decays very rapidly with the length, a simulation for a 
large system may still be confined to the sector of zero 
winding number. 



B. Spin and dimer correlations 

In order to characterize the different phases realized by 
the CAP states, we evaluate order parameters for detect- 
ing antiferromagnetic (Neel) order and VBS order. Neel 
order can be characterized using the standard two-spin 
correlation function, 



C(ry) = <Sr j -S Pj )(-l)^«+»«) 1 



(15) 



where we use ry to denote the vector separating the lat- 
tice sites i and j and the phase factor cancels the signs of 
the staggered spin correlations obtaining in the systems 
we study. Alternatively, one can study the full sublattice 
magnetization averaged over the whole system; 



(16) 



where 4 = +1 on sublattice A and <f>i = —1 on sublat- 
tice B. Since the singlet AP and CAP states manifestly 
cannot break the spin-rotation symmetry, order must be 
detected in the squared order parameter, (m^), which 
in the limit of large system size will be identical to the 
long-distance spin correlation (|15[) . 

To accurately locate an antiferromagnetic phase tran- 
sition, the Binder cumulant is very useful. It is defined 
according tc42 



2 V 5 (m 



(17) 



where the factors are chosen for the 3-component Neel 
order parameter such that U(N oo) = in the disor- 
dered phase (where the order-parameter distribution is a 



Gaussian with zero average) and U(N — s- oo) = 1 in the 
ordered phase (where the radial distribution is peaked 
at non-zero m s ). Typically crossing points of U graphed 
versus a control parameter for different system sizes ap- 
proach the critical point vary rapidly as a function of 
increasing system size. 

To characterize VBS order we use the dimer correlation 
function, defined as 



D xx {vij) = {B x {Ti)B x (Tj)), 
D yy( r ij) = ( B y(n)B y (r 3 )), 



in terms of the bond operators 

B x {vi) = S ri ■ 
B y (r t ) = S ri ■ 



'l-j+XJ 



'r;+y i 



(18) 



(19) 



directed along the unit lattice vectors x and y. We will 
not need the mixed x-y correlations here. In some cases 
we will characterize VBS order by the long-distance be- 
havior of (|18|) . The states we will be studying have a 
2-site VBS unit cell, forming a staggered weak-strong- 
weak-strong pattern in one dimension and an analogous 
columnar pattern in two dimensions. In both cases we 
can extract the dominant component of the correlations, 
corresponding to the squared order parameter, by taking 
the appropriate difference of (JTHJ) evaluated at nearby 
distances. We here use a symmetric version of this dif- 
ference; 

D* x (t) = D xx (r) - I [D xx (r - x) + D xx (r + x)] , (20) 

and a function D* y (r) for y-oriented dimers defined anal- 
ogously. We will also study the full order parame- 
ter, which in two dimensions can be defined using the 
q = (it, 0) and q = (0, tt) Fourier transforms of the 
nearest- neighbor bond correlations (1191) : 



i 



(21) 



The magnitude D of the order parameter can be com- 
puted as the square-root of the average squared opera- 
tor, (D 2 ) = (D x ) + (D 2 ). In addition to the expectation 
values, we will also investigate the probability distribu- 
tion P(d x ,d y ), in which emergent U(l) symmetry can 
be detected. Here d x and d y are the expectation values 



of the corresponding operators (|2ip evaluated in a given 
sampled configuration based on the transition graph. We 
refer to Ref. [52j for further details on this quantity, which 
is not a conventional quantum mechanical expectation 
value but still very useful for characterising VBS states. 

All of the above two- and four-spin correlations are 
related to the transition-graph loops generated in the VB 
MC sampling process. For instance, the estimator for the 
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two-spin correlation is given by& 



3.4 



(V a \S ri -S r M 



(V a \V p 



=c 4' 

0. 



(22) 



where and denote sites i and j belonging to the 
same loop and different loops, respectively, and the sign 
in the case is + and — for spins on the same and 
different sublattices, respectively. From Eq. (|22|) one can 
also obtain a very simple expression for the estimator for 
the squared staggered magnetization, 



(v a \v ) 



(23) 



where Ci is the size (the number of sites) of loop £ . 

Both the dimer correlation function and the fourth 
power of the staggered magnetization involve four-spin 
correlations. Detailed descriptions on how to calculate 
these based on the transition graph of two VB configura- 
tions can be found in Refs. and [U Here we only write 
down the expression for the fourth power of the staggered 
magnetization, needed for the Binder cumulant (|17[) . 



(Va\V ) 



E^(E^) -sE4 (24) 
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which is also solely determined by the sizes of all loops 
formed in the transition graph. We note that the Binder 
cumulant of the VBS order parameter is much more diffi- 
cult to evaluate, since its definition in analogy with (|17p 
requires eight-spin correlations. While these also in prin- 
ciple can be evaluated in terms of the transition-graph 
loops r the expressions are quite complicated to imple- 
ment in practice and we have not done so. 
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FIG. 3: (Color online) The upper panel shows the crossing 
behavior of the Binder cumulant U(k) defined in Eq. (|17|) 
for several different chain lengths L when A = 1. The ap- 
proach to 1 for small k and for large k corresponds to the 
presence and absence of Neel order, respectively. The cross- 
ing point approaches the critical value of k. The lower panel 
demonstrates extrapolations to the thermodynamic limit of 
the critical k c by fitting crossing points of (L, 2L) pairs to the 
power-law correction 



III. NEEL TO VBS TRANSITION IN ONE 
DIMENSION 

In one dimension, the standard AP states given in Eq. 
(|6]) are able to reproduce a Neel- VBS transition without 
correlation factors. We will study this ID transition care- 
fully in this section, using the very efficient loop update 
of the VB configurations. 

It is natural to study the evolution of the state as a 
function of some parameter governing the long-distance 
behavior of the amplitudes, e.g., using the power law 
h(r) = l/r K with tunable n or an exponential form. Here 
we will use the power-law. However, it is already known 
that the nature of the state is not just determined by 
the asymptotic behavior of h(r), but also depends on 
details of the short-bond weightsJ^ In addition to the 
exponent k we here tune the shortest-bond amplitude 
h(r = 1) = A. The wave function is, thus, explicitly 
given by 

/ i \n r (V) 

#)^ ni(v) nb . ( 25 ) 



where n r (V) again refers to the number of bonds of 
length r in the bond configuration V. 

It is clear that for A > and large k this AP state 
is a VBS, since in the limit n — > oo only two configura- 
tions contribute; those with r = 1 bonds on alternating 
links. For small k there is instead Neel order but no VBS 
order Note that long-range order corresponding to bro- 
ken SU(2) symmetry is possible in this kind of ID system, 
since viewed as a classical statistical-mechanics problem 
there are long-range interactions (since the bonds have 
unbounded length), and the Mermin- Wagner theorem^ 
prohibiting ID Neel order does not apply. Note also again 
that the AP wave function is a singlet and, thus, the 
SU(2) symmetry is not actually broken (as in any calcu- 
lation targeting the singlet ground state) . The magnitude 
of the Neel order measured by (m^), Eq. (|16[) . or the long- 
distance correlation function (fT5j) can still evolve toward 
a non-zero value as the system size grows, tending to the 
square of the symmetry-broken value of m s in the corre- 
sponding thcrmodynamic-limit state with no constraint 
on the total spin. 
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FIG. 4: (Color online) Phase diagram of ID AP states with 
tuning parameters k and A, as defined in Eq. (|25fl . The cir- 
cles are calculated transition points and the curve is a guide 
to the eye representing approximately the boundary between 
the long-range ordered Neel (below) and VBS (above) phases. 
The inset exemplifies long-distance spin correlation functions 
inside the phases and at the critical point when A = 1; the 
black squares correspond to k = 1.6 (inside VBS phase) and 
the green triangles are for k = 1.4 (in the Neel phase). The 
red circles show the behavior at the critical point. 

Beach has previously studied Neel ordering in this 
class of wave functions (with a somewhat different 
paramctrization of the short-bond amplitudes)^ He 
found a continuous transition between the Neel state and 
the non-magnetic state. Here we also investigate the VBS 
correlations and find a single transition point where both 
the spin and dimer correlations are critical. We study the 
evolution of the transition in the plane (k, A). 

For fixed A, in order to find the critical value of k c of 
the AP state we study the Neel Binder cumulant ([TT]) . 
The behavior of curves for different system sizes L cross- 
ing each other as a function of k is illustrated in the upper 
panel of Fig. |3l The crossing points do not fall exactly 
on a single point due to subleading size corrections. We 
observe a systematic smooth drift of the crossing points 
as the system size is increased. In order to eliminate this 
size effect and determine the critical point from data such 
as those in Fig. [31 we extract K-valucs corresponding to 
crossing points of (L, 2L) size pairs, and plot these points 
against 1/L, as shown in the lower panel of Fig. [3] Wc 
then extrapolate these values to L — > oo and obtain hl c . 
The fitting function we use here for extrapolation is the 
standard power-law^ 

f c (L,2L) = Kc + -^. (26) 

The extrapolated k c values versus A are plotted in Fig. |U 
the phase diagram of ID AP states with the two tuning 
parameters A and n. The inset of this figure demonstrates 
the qualitatively different behaviors of the spin correla- 
tion functions in the two phases and at the critical point, 
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FIG. 5: (Color online) Staggered spin-spin (upper panel) and 
dimer-dimer (lower panel) correlations of ID AP states at the 
largest distance, graphed versus the chain length at c = k c 
for different short-bond amplitudes A. All lines are fits to the 
form aL~ b . 



using A = 1 results as an example. At k = 1.6 the 
correlations decay faster than power-law, as is expected 
for a non-magnetic VBS ordered state. In contrast, at 
k = 1.4, the correlations for small L first decay some- 
what but then converge to a non-zero value for larger 
L (even increasing somewhat for large systems), demon- 
strating the presence of long-range Neel order. At the 
critical value k c extracted using the Binder crossings as 
explained above, the decay of the correlations arc consis- 
tent with a critical, power-law form. 

To determine whether the VBS correlations are also 
critical at the k c points extracted from the Neel Binder 
cumulant, we further study both the spin and dimer cor- 
relations at these points. The results confirm the expec- 
tation of a common critical Neel and VBS point. By 
studying chains as large as L = 4096, we can extract 
the exponents governing the critical correlation functions 
with relatively small error bars (thanks to the powerful 
VB MC loop update discussed in Sec. |H]). The analysis 
of the power laws is presented in Fig. [5] Note that in 
order to avoid boundary modifications of the power-law 
correlation functions as a function of the distance r in 
systems of fixed L, we study the long-distance correla- 
tions versus the system size, with r = L/2 for the spin 
correlations and the staggered component of the dimer 
correlations extracted based on r = L/2 and L/2 — 1 
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FIG. 6: (Color online) The continuously varying spin (a) and 
dimer (/3) decay exponents of the ID AP state (J25J as a func- 
tion of the short-bond amplitude A. The exponents corre- 
spond to the power-law decay of the correlation functions; 
C(r) ~ r~ a , D*{r) ~ r~ p . The points are calculated values 
and the curves are guides to the eye. 



data according to Eq. (|2"U)) [where it should be noted that 
D(L/2 - 1) = D(L/2 + 1) for a periodic chain]. 

As A increases, larger system sizes are needed to ob- 
serve the asymptotic critical forms. Especially for the 
largest A studied, A = 8, one can observe in Fig. [5] (up- 
per panel) a clear cross-over of the spin correlation func- 
tion from a rapidly decaying short-distance form to the 
asymptotic power-law form. The straight lines in Fig. [5] 
are fits to the simple asymptotic form aL~ b . We have 
also tried to include shorter chains in an analysis includ- 
ing corrections, by fitting to the form aL b + cL d . This 
form is, however, not capable of describing the small size 
effect in this model (in contrast to 2D critical spin liquid 
RVB states, where this form works very well^i). In any 
case, the large-L behaviors appear to be reasonably well 
converged to the simple power law and the exponents ex- 
tracted should be reliable. An exception is A = 0, for 
which the dimer correlations decay very rapidly and are 
too noisy to allow the exponent (3 to be reliably deter- 
mined (and we have therefore not graphed these correla- 
tions in Fig. [5]). It is even possible that the VBS state 
for A = is of a different kind than for A > 0. Further 
studies will be needed to settle this issue. 

We plot the extracted critical exponents as a function 
of A in Fig. [6] The exponents vary continuously with A, 
with the dimer exponent decreasing monotonically and 
the spin exponent increasing. An interesting conclusion 
that can be drawn from these results is that the criti- 
cal state becomes increasingly "quasi VBS ordered" with 
increasing A, with the decay exponent of the dimer cor- 
relations perhaps vanishing as A — > oo, although this 
is difficult to confirm definitely (because the simulations 
become increasingly difficult for large A). The behavior 
is in line with the expectation that a large lambda fa- 
vors VBS ordering because of the predominance of the 
very shortest bonds, i.e., when moving on the critical 



line toward higher A the density of short bonds increases, 
and this leads to a strengthening of the VBS quasi-order. 
At the same time, the exponent of the spin correlations 
appear to approach I. However, Neel order still exists 
for large A when reducing n from the critical value. In 
terms of the transition graph estimators of the correlation 
functions, VBS correlations correspond to certain loop 
correlations^ while Neel order is related to the presence 
of long (cx L) loops. While long-range Neel and VBS or- 
ders are mutually exclusive in these states, the Neel state 
in the neighborhood of the critical curve for large A ap- 
proaches a coexistence situation. Here the magnitude of 
the Neel order parameter also becomes very small, how- 
ever, and the coexistence is therefore not robust. 

The VB formulation of the ground state can be viewed 
as a ID classical statistical mechanics problem, but at the 
same time it should also correspond to a path-integral 
formulation in 1 + 1 dimensions (with some underlying 
parent Hamiltonian) . One may then expect the system 
to be classifiable according to the standard 2D conformal 
field theories by a central charge c. Varying critical expo- 
nents, as we have found here, normally imply c > 1, but 
the fact that the system includes long-range interactions 
may invalidate this requirement, although it is not clear 
how the power-law bond length translates into effective 
interactions in an underlying parent Hamiltonian (and 
the interactions in it may well be short-ranged). One no- 
table aspect of the AP states is that they are not able 
to reproduce the ground state of the critical Heisenberg 
chain, where a = j3 = We will address this issue 
further in Sec. [V] with variational AP calculations for the 
Heisenberg chain. 

It would be interesting in the future to compute the 
bipartite entanglement entropy of the ID AP states, to 
test its system size scaling and consistency between c 
extracted from i t 45 ' 46 and from the correlation functions. 
Such calculations can also be carried out using the VB 
MC sampling scheme used here i 23 ' 47 



IV. NEEL TO VBS TRANSITION IN TWO 
DIMENSIONS 

The Neel state is known to be the ground state of the 
square lattice Heisenberg antifcrromagnct with homoge- 
neous nearest-neighbor couplings. There is convincing 
numerical evidence^ as well as mean-field arguments 1 ^ 
showing that the standard AP states with power-law de- 
caying amplitudes of the asymptotic form h(r) = 1/r 
is an optimal variational wave function for the 2D Neel 
state (for any finite size, where the ground state is a sin- 
glet with no explicitly broke symmetry). 

The AP states have no Neel order for rapidly decay- 
ing (exponentially or according to a power law 1 /r K with 
large k) bonds.— An extreme case is the state that con- 
tains only nearest neighbor bonds (dimers). Such a short- 
bond VB state on the square lattice normally corresponds 
to an U(l) spin liquid with critical VBS correlations^ 1 ^ 
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FIG. 7: (Color online) Results for the CAP states with the 
case (I) parametrization. (a) The Binder cumulant of the 
staggered magnetization as a function of the weight p favor- 
ing parallel dimer alignment (ci = p in Fig. [JJ with the other 
configurations suppressed by setting Ci>i = 1/p)- The grow- 
ing negative peaks of the cumulant indicate the location of the 
first-order phase transition developing as a function of the sys- 
tem size, (b) The columnar component D*(r) of the dimer 
correlation function at the largest distance, r = (L/2, L/2), 
plotted against p. The correlation function approaches zero 
for large systems in the non-VBS state and becomes finite 
in the VBS. The behavior of the correlation function tending 
to a step-function as L increases, and the curves for differ- 
ent L crossing each other, is in accord with the first-order 
phase transition signaled by the pronounced negative cumu- 
lant peaks in (a). 



in contrast to the ID AP VBS state discussed in the pre- 
vious section. One should expect the 2D AP state to 
turn into a long-range ordered VBS when appropriate 
bond correlations are included. On the square lattice, 
which we will consider here, the simplest kinds of VBS 
states (with 2-site unit cell) form columnar, staggered, 
or plaquette ordering patterns. 

We here study the Neel to VBS transition on a square 
lattice within the CAP states by imposing bond correla- 
tions that favor or suppress only certain types of short- 
dimer alignments. All possible configurations of short 
dimers connected to a pair of nearest-neighbor sites on 
a square lattice are shown in Fig. [TJ We assign a weight 
to each of those two-dimer configurations according to 
Eq. fTJ), with all weights with ri,r2 ^ 1 set to 1. To 
simplify the notation we here use c,; for the short-bond 
correlation factors, instead of C(ri,r2), with the corre- 
spondence between the two shown in Fig. [TJ For the 
special case of there being a single bond connecting the 
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FIG. 8: (Color online) Finite size scaling behavior of the 
Binder cumulant minimum of the CAP states, case (I). The 
negative minimum increases linearly with L 2 , indicating a 
first-order transition in this case. 



two reference sites we set Co = 1 as a normalization factor 
for the correlations. 

To reduce the number of control parameters in our 
simulations we introduce a single parameter p such that 
Ci = p > 1 for favored two-dimer configurations i, while 
Ci = 1/p < 1 for unfavored configurations and p = 1 for 
cases that are considered "neutral" . If all dimers are un- 
corrected, i.e., p = 1, the state reduces to the standard 
AP state, for which we choose the single-bond amplitudes 
to be h(r) = 1/r 3 . This choice, which we keep also for 
the CAP states, is motivated by the fact that this gives 
the correct description of the Neel state. The generalized 
CAP states we use in our simulations are therefore char- 
acterized by the single parameter p controlling the bond 
correlations. 

Below we investigate two different parametrizations of 
the bond correlations. In both cases we use c\ = p > 1 
to locally favor the columnar or plaquette VBS pattern 
(and whichever of these two VBS patterns that actually 
will be realized is not clear from the outset). Other types 
of dimer correlations are suppressed, by setting C2 = C3 = 
C4 = 1/p in the first case — case (I) — while they are set to 
neutral, C2 = C3 = C4 = 1, in case (II). 

We destabilize Neel order by increasing the control pa- 
rameter p and study the phase transition into a VBS. 
For case (I), we have found a first-order Neel to colum- 
nar VBS transition, while for case (II) we have found a 
continuous transition into a critical U(l) spin liquid, fol- 
lowed by a second continuous transition to the columnar 
VBS. We discuss the two cases in order. 



A. A first-order Neel to VBS transition 

Case (I) again corresponds to favoring parallel VB 
bond configuration by setting c\ = p > 1 in Fig. [TJ and 
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suppressing fluctuations by setting Cj>i = 1/p. As in the 
ID case discussed in Sec. lIIIi we here first use the Binder 
cumulant of the staggered magnetization U to detect the 
transition of the Neel order. This quantity is also useful 
for distinguishing between a first-order and continuous 
phase transitions. As shown in Fig. Eta), the Binder cu- 
mulant as a function of the control parameter p exhibits 
a minimum separating the Neel phase, where U — > 1, 
and a non-magnetic phase, where U —¥ 0. The minimum 
value of U is negative for all system sizes we studied, and 
the negative peak becomes narrower and deeper as the 
system size increases. In fact, the negative peak diverges 
as — L 2 when L — > oo, as plotted in Fig. [HI which provides 
strong evidence^ for a first-order phase transition. 

Note that the divergence of U m i n of the form L d ex- 
pected for a classical d-dimcnsional system could in prin- 
ciple change to L d+Z for a quantum system, where z 
is a dynamical exponent^ However, our definition of 
the Binder cumulant is purely a real-space definition and 
does not include integration over the imaginary time di- 
mension (which we do not even have access to because it 
relics on a path-integral formulated using the underlying, 
unknown parent Hamiltonian). The form L 2 therefore is 
expected. 

To check whether the non-magnetic phase exhibits 
VBS order, we next compute the columnar component of 
the dimer correlation defined according to (|20|) . Fig. [7th) 
shows D*(r) for the largest distance, r = (L/2, L/2). 
The correlation indeed converges to a non-zero value in 
the non-magnetic phase, tending to a step function at p c 
as L increases. The location of the discontinuity coin- 
cides with the point where the Binder cumulant reaches 
its minimum in Fig. [T^a) . Note also that the curves for 
different L cross each other. This size-independence of 
the order parameter (as supposed to size-independence 
after multiplying with some power of L corresponding to 
an exponent of critical correlations) at the transition also 
supports a first-order scenario. 

The location of the transition point p c in the thermo- 
dynamic limit can be determined, e.g., by extrapolating 
the U minimum location p m i n (Fig. [7|) to the infinite-!/ 
limit. The finite-size scaling plot in Fig. [H^a) shows that 
the finite-size shift of the transition point defined in this 
way is consistent with oc L~ 2 for large L, where the ex- 
ponent 2 again is the one expected based on scaling at 
a first-order transition, as discussed above. We estimate 
p c « 1.500 from an extrapolation to L — > oo. 

In the regime for p < p c the cumulant for different sys- 
tem sizes exhibits crossing points versus p. We expect 
that the crossing points should coincide with the mini- 
mum location when L = oo. By finite-size extrapolation 
of crossing points p* for pairs of two system sizes L/2 
and L, shown in Fig. G^b), we estimate p* « 1.500 in 
the thermodynamic limit, in perfect agreement with the 
result obtained from the cumulant minimum. 

In addition to the Binder cumulant signaling the tran- 
sition of the Neel order, we also estimate the transition 
point of the VBS order from the scaling of the crossing 
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FIG. 9: (Color online) Finite-size scaling plots for extracting 
the location of the phase transition in CAP case (I), (a) the 
location p m i n at which the Binder cumulant reaches its min- 
imum, (b) the crossing point of the Binder cumulant for sys- 
tem sizes {L/2, L), and (c) the crossing point of the columnar 
dimer correlation function. From finite-size extrapolations of 
these three quantities, assuming 1/L 2 dependence, we obtain 
consistently the transition point p c located within the range 
1.500 — 1.505. Note that there are still some visible deviations 
from the assumed form and a more precise determination of 
p c would require data for still larger systems. 



points p' of the long-distance dimer correlation function. 
The inset of Fig.[7tb) shows a zoom- in of the region where 
the crossings occur. As shown in Fig.[SJc), these crossing 
points also appears to shift as p' ~ L~~ 2 for large L, and 
the extrapolation to L — > oo yields an estimated location 
of the transition point p' c ps 1.503. This is marginally 
above the two other estimates discussed above, but given 
the very small range of data points for which the L~ 2 fits 
work well (we have used the points for the three largest 
systems in all cases, but the data still show some non- 
asymptotic curvature here), this result is still consistent 
with a single Ncel-VBS transition point. Larger system 
sizes would be required to extract the location of this 
point more precisely. 

Finally for case (I) we examine the histogram P(d x ,dd) 
of the order parameters d x and d y , corresponding to VBS 
order with x- and y-oricntcd dimcrs, Eq. (|21[) . In Fig. 1101 
we show P(d x ,d y ) for L = 28 at p = 1.46, 1.48 and 1.50. 
At p — 1.46 inside the Neel phase, the distribution has 
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FIG. 10: (Color online) VBS order parameter distribution 
P(dx,d y ) for L = 28 CAP states in case (I). Brighter re- 
gions correspond to higher density. Left panel: The distri- 
bution at p — 1.46 is circular-shaped with a central peak, 
showing U(l) symmetry in the Neel state. Right panel: At 
p — 1.50 there are four peaks at the Z4-symmetric angles 
cj> = aicta,n(D x /D y ) = 0, 7r/2, tt, 3tt/2, reflecting columnar 
VBS order. Middle panel: At p = 1.48, in the transition 
region, the 5-peak distribution shows coexistence of the Neel 
and VBS order, providing evidence for a first-order transition. 



a circular shape with a central peak. At p = 1.48 in 
the transition region, the distribution shows coexistence 
of the Neel order (characterized by the central circular 
region) and the columnar VBS order [characterized by 
the four narrow peaks at angles (f> = arctan(d y /cic) = 
0, 7r/2, tt, 3tt/2]; this again provide clear evidence for a 
first-order Neel- VBS phase transition. At p = 1.50, only 
barely inside the columnar VBS phase of this finite sys- 
tem, the distribution exhibits only the four VBS peaks. 

From the many consistent results discussed in this sub- 
section, we can conclude that the CAP states with fa- 
vored parallel dimer pairs and suppressed flips of such 
pairs can characterize a first-order phase transition be- 
tween the Neel and the columnar VBS phases. In such a 
CAP state, the Neel order is destroyed by formation of 
parallel dimcrs as the weight p increases. A first-order 
transition is also of course what would normally be ex- 
pected for an order-order transition involving two unre- 
lated order parameters. Note also that the system sizes 
studied in this section were rather modest; up to L = 36 
(while much larger systems, up to L = 128, will be con- 
sidered in the next section). For larger systems it be- 
comes very difficult to obtain good statistics and smooth 
curves versus the control parameter, because of hysteresis 
effects related to the first-order nature of the transition. 
Still, as we have shown, the system sizes studied are suf- 
ficient to study the asymptotic finite-size scalimng forms. 

In the context of VBS ground states of Hamiltonians, 
a first-order transition was previously observed with a 
J-Q model with the multi-spin interaction Q arranged 
to favor a staggered stated In that case, the first-order 
transition was expected, because local dimer fluctuations 
are strongly suppressed with this kind of bond order. Al- 
ternatively, one can make an argument based on the na- 
ture of vortex-like defects in the VBS^i In contrast, in 
a columnar state parallel short-dimer pairs can fluctuate 
by 90° rotation, unless such fluctuations are energetically 
expensive. In the DQC theory, these fluctuations are es- 
sential and correspond to an emergent U(l) symmetry of 




FIG. 11: (Color online) Results for CAP states, case (II), for 
different system sizes plotted versus the control parameter 
p. (a) The binder cumulant of the staggered magnetization 
exhibits a shallow peak at the transition from the Neel to 
a non-magnetic state, with p c \ ~ 1.28. (b) The columnar 
dimer correlation function indicates a second transition, into 
the VBS, at p a i « 1.65, with a completely disordered phase 
for Pd < p < p c l2. 



the VBS order parameter, which has been confirmed in 
J-Q models with plaquette VBS ground states j^^^^ 
In the CAP state considered here, we suppressed the fluc- 
tuations out of the perfect columnar state by the use of 
correlation weights and the observed first-order transi- 
tion then is in line with expectations based on the DQC 
theory and the earlier studies of the Neel- VBS transi- 
tions in various models. 



B. A critical U(l) spin liquid and a second-order 
transition to the VBS 

Given the findings and discussion in the previous sec- 
tion, we now investigate whether the removal of the corre- 
lation factors suppressing dimer fluctuations can change 
the nature of the Neel- VBS transition of the CAP states. 
In this case (II) we only set c\ = p > 1 and keep the other 
correlation factors in Fig. [T] as neutral; c% = ca = c& = 1. 

Fig. [TTJ shows the Neel Binder cumulant and the stag- 
gered dimer correlation function against the control pa- 
rameter p. Like in case (I), the Neel order, characterized 
by the Binder cumulant tending to 1 as L grows, survives 
in the small p region up to a certain value of p, and a sub- 
stantial columnar dimer correlation D* sets in when the 
Neel order is destroyed. 
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FIG. 12: (Color online) Finite size scaling behavior of the 
Binder cumulant minimum in case (II). The negative mini- 
mum (obtained by interpolation of the data shown in Fig. 1 1 1 j) 
increases only logarithmically with the system size L, indicat- 
ing a continuous transition. 



We notice that the negative peak of the Binder cu- 
mulant, which occurs only for large systems, is less pro- 
nounced than the cumulant peak in case (I). A nega- 
tive Binder cumulant is often taken as evidence for a 
first-order transition^ but there are now known exam- 
ples of rigorously understood continuous classical phase 
transitions associated with this behavior, e.g., the 2D 4- 
state Potts model and the related Ashkin- Teller model 
in the neighborhood of this special point^ In such cases 
of "pseudo-first-order" scaling the minimum U m i n of the 
Binder cumulant diverges much slower than the expected 
L d form at a first-order transition (or, in some cases, pos- 
sibly even converges to a finite value) . 

The finite-size scaling plot Fig. [T2l shows that — U m i n 
for CAP in case (II) grows only logarithmically with 
the system size. Thus, we conclude that the Neel or- 
der here vanishes in a continuous transition associated 
with pscudo-first-order behavior (related to anomalies in 
the critical order-parameter distribution). 

Next we determine the phase boundaries more pre- 
cisely. For a continuous phase transition, a frequently 
used method to determine the critical point is to find the 
unique asymptotic crossing point of the Binder cumulant 
for different system sizes. Due to finite size effects, the 
cumulant curves will intersect at a single point only when 
the system sizes are sufficiently large. For our case, we 
find the intersection point of the cumulant for L > 80, 
and it is located at p ss 1.276 [Fig. flUT a)]. which can be 
identified as the Neel-spin liquid transition point p ic . As 
the cumulant for a large system size exhibits a negative 
peak before it vanishes in the spin liquid phase, we also 
extrapolate the location of the peak to L — > oo to deter- 
mine the critical point from another route. By doing so, 
we estimate p\ c « 1.276, in perfect agreement with the 
cumulant intersection point. 

Beyond the order of the transition, another major dif- 
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FIG. 13: (Color online) (a) The Binder cumulant for large 
system sizes in case (II). The curves intersect at a point, p w 
1.276, which separates the Neel phase and the spin liquid, (b) 
Finite size scaling of the location p m i n at which the cumulant 
reaches its minimum. The critical point of the Neel-spin liquid 
transition in the thermodynamics limit is also here estimated 
as pic ss 1.276, from an extrapolation to L — oo 



ference from case (I) is the behavior of the dimer corre- 
lations. As seen in Fig. UTT b). there is a wide interme- 
diate region between the Neel phase and the VBS phase 
for larger p (where the correlation clearly converges to a 
non-zero constant for large systems); in this intermediate 
region, the dimer order parameter still decays versus the 
system size and, for the largest systems, the curve ver- 
sus p develops a behavior suggestive of a second phase 
transition between a disordered state and the VBS above 
p = 1.6. The dimer correlations decay in the intermediate 
region as a power law, D*(r) ~ r - ^, with the exponent 
j3 depending on p, as shown in Fig. [TJ] Algebraically de- 
caying dimer correlations were previously found in short- 
bond resonatin g va lence bond (RVB) spin liquids inves- 
tigated in Refs. I2lll22l . We thus tentatively identify this 
intermediate p-region with critical dimer correlations as 
a spin liquid in the same class of RVB states, for which 
it is known that the exponent of the dimer correlations 
depends on details of the bond fugacitics^i 

We next study the distribution of the VBS order pa- 
rameter P(d Xl d y ). The examples of distributions shown 
in Fig. [15] for L = 64 are ring shaped in the intermedi- 
ate region (at p = 1.4 and p = 1.5) before evolving into 
the expected Z4-symmetric shape in the larger-p regime 
where columnar VBS order is formed (as seen in the loca- 
tion of the peaks; a 45° rotated distribution would corre- 
spond to a plaqucttc VBS). The same kind of ring-shaped 
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FIG. 14: (Color online) Long-distance columnar dimer corre- 
lation function versus system size for different p in the inter- 
mediate region, the correlations decay as a power-laws (with 
the decay exponent /3 = 0.57 for p = 1.4, 0.43 for p — 1.5, and 
0.32 for p = 1.6), in the spin liquid phase and are size inde- 
pendent in the VBS state (p = 1.8). At p — 1.7 the behavior 
is approximately a power-law but with a slight up-turn for 
the largest size, indicating a weakly VBS ordered state here. 

P{d x ,d y ) distribution was also found in the prototypi- 
cal short-bond RVB spin liquid (although only when the 
bond configurations are restricted to the dominant topo- 
logical sector of zero winding number) i 21 i 22 

An important issue here is whether the VBS hosts 
emergent U(l) symmetry when the critical point is ap- 
proached. To investigate this, we also need to determine 
the location of the point where the VBS becomes long- 
range ordered. This is not easy to do just based on the 
scaling behavior versus the system size in Fig. [14l because 
there is a whole critical VBS phase and the change from 
the power-law decay to convergence to a samall non-zero 
constant is subtle. We can instead characterize the VBS 
state by the quantity cos(40), where <j) = arctan(c?j, / d x ) is 
the angle in the VBS order parameter (d x , d y ) computed 
for the individual VB configurations, i.e., based on the 
histogram, P(d x ,d y ). This expectation value measures 
the degree of the developed Z4 symmetry of the VBS or- 
der parameter. In the spin liquid (as well as in the Neel 
state), where the U(l) symmetry is preserved as L — >• 00, 
we have (cos(4(/>)) = 0. In the VBS states in which the 
distribution for large systems develops Z4 symmetry, we 
have (cos(4</>)) — > 1 as L — > 00 for a columnar VBS (while 
it approaches —1 for a plaquette VBS). As cos(4^>) is 
dimensionless, one expects (cos(4</>))-curves for different 
system sizes to asymptotically become size-independent 
at the transition point into the VBS state. 

As seen in Fig. [ToTa), for the case- (II) CAP, a cross- 
ing point develops at a non-zero value of (cos(40)), at 
p m 1.65, which we thus identify as the liquid- VBS tran- 
sition point. The fact that cos(4</>) > at this point 
shows that there is no emergent U(l) symmetry at the 
VBS-liquid transition, since the order parameter remains 
Z4 symmetric exactly at the transition point. For ref- 



FIG. 15: (Color online) Histograms of the VBS order param- 
eter defined in Eq. (|2ip shown for L — 64 systems in CAP 
case (II). In the spin liquid phase with p = 1.4 and p — 1.5 
[(a) and (b), respectively], the distributions are ring-shaped 
with weight at all angles (bright regions). As p increases to 
p = 1.6 (c) and p — 1.7 (d), the U(l) symmetric distributions 
evolve into Z4-symmetric ones, with higher densities at the 
angles 0, 7r/2, it, 3n/2 corresponding to columnar VBS order. 



crence, in Fig. [T7] we show QMC results for the J-Q 
model with 6-spin columnar Q-interactions, for which the 
order-parameter symmetry was previously analyzed in a 
slightly different way^ The results for (cos(4<^)) here 
show crossing points decaying toward in the vertical 
direction. The system sizes are not yet sufficiently large 
to see that the crossings tend toward the critical point, 
(J/Q) c ~ 0.66. The contrast with the CAP states in 
Fig. lTBT a) is stark, however, with the absence of four-fold 
symmetry — presence of emergent U(l) symmetry — at the 
transition being very plausible. 

Since the results shown in Fig. 116( a) appear to give 
a rather precise estimate for the transition point, p2 C = 
1.650(5), without any scaling needed of the vertical axis, 
we now use this result to investigate scaling of the VBS 
order parameter. In Fig. 1181 in order to achieve data 
collapse, we have rescaled the long-distance dimer corre- 
lation functions in Fig. I13f b) both vertically, multiply- 
ing by L@ , with (3 = 0.22, and horizontally, multiplying 
(p — P2c) by L 1/v , with \jv = 0.44. The correlation- 
length exponent is, thus, v sa 2.3, which is anomalously 
large and definitely rules out a first-order transition (in 
which case v = l/d= 1/2 would be expected). 

The (cos(40)) curves in Fig. [TBTa) can also be scaled 
in the horizontal direction to collapse all the data onto 
a single curve, as shown in panel (b). In cases where 
there is emergent U(l) symmetry, this procedure, using 
the control parameter scaled as L x l av (j> —p c ), gives the 
correlation-length exponent v multiplied by a number 
a > 1,— reflecting the faster divergence of the length- 
scale A controlling the emergent symmetry; A ~ £ a . In 
the case at hand, we have already concluded that there 
is no emergent U(l) symmetry, since (cos (40)) remains 
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FIG. 16: (Color online) The quantity (cos(4<?!>)) with </> = 
arctan(d y / d x ) measures the degree of Z4 symmetry in the 
VBS order parameter and also gives a good estimate of 
the location of the liquid- VBS transition. In the spin liq- 
uid the distribution P(d x ,d y ) is U(l) symmetric, leading to 
(cos(40)) = for L — > 00, while in the columnar VBS phase, 
where the distribution is Z4-symmetric, (cos(40)) approaches 
1. (a) The curves for different system sizes cross at a single 
point located at p = 1.65, which is identified as the transition 
point; P2c = 1.650(5). Since (cos(4</>)) > at this point, there 
is no emergent U(l) symmetry at this transition, although it 
appears to be a continuous transition based on other results, 
(b) When p — p2c is scaled with the system size L raised to the 
power L 1 / au , with av « 1.82, the curves for different systems 
collapse onto each other. 



non-zero as p — > p2c, and the exponent av, with a < 1 
should instead reflect a shorter length-scale (irrelevant 
operator), governing the reduction of the angular VBS 
fluctuations ((cos(4t/>)) — > 1) in the VBS and the growth 
of these fluctuations ((cos (40)) — > 0) in the liquid. We 
find very good scaling, and indeed the factor a ~ 0.8 is 
clearly less than one. 

We finally investigate the nature of the spin correla- 
tions in the spin liquid phase. Results for the squared 
sublatticc magnetization for several points representing 
the three different phases arc shown in Fig. [19] Here 
we plot the results on a log-log scale, in order to study 
power-law correlations. In the Neel state the sublatticc 
magnetization approaches a non-zero constant, while in 
the liquid and VBS states we observe a clear 1/L 2 decay. 
This is the form expected with exponentially decaying 
spin-spin correlation functions. We see a power-law be- 
havior with a non-trivial exponent, ~ L~ a , with a = 1.55 
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FIG. 17: (Color online) The degree of Z4, anisotropy of the 
VBS order parameter of the J-Q model with 6-spin inter- 
actions for different system sizes. Here the curve crossings 
tending to (cos (40)) = with increasing L demonstrates the 
emergent U(l) symmetry at the Neel-VBS transition of this 
model. 
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FIG. 18: (Color online) The data of Fig HUb) in the neigh- 
borhood of the liquid- VBS transition, rescaled to extract the 
exponent ft of the dimer correlation function, D*(r) ~ r~ fl , 
here with ft — 0.22, and the correlation length exponent v, 
here with v « 1/0.44 w 2.27. 



only at the Neel-liquid critical point. These results again 
confirm that the liquid is of the same type as the proto- 
typical AP RVB states (i.e., with no correlation factors), 
where the varying power-law for the critical VBS corre- 
lations corresponds to a varying stiffness constant in a 
mapping to a height modc h 21 i 55 



V. VARIATIONAL CALCULATIONS 

In this section we explore variational optimization of 
CAP states, carrying out energy minimization based on 
derivatives along the same lines as in Refs. [l9H20l . We con- 
sider two models. First, the standard Hcisenberg model, 
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FIG. 19: (Color online) Finite-size dependence of the squared 
staggered magnetization for different values of p. In the Neel- 
phase (p — 1.260), m 2 converges to a non-zero value for L — > 
oo. At the Neel-spin liquid critical point, p = 1.276, it scales 
as m 2 ~ L -1 ' 55 . We see m 2 ~ L~ 2 inside the spin liquid phase 
(p — 1.4), at the liquid- VBS critical point (p = 1.65), as well 
as in the VBS phase (p — 2), indicating an exponentially 
decaying spin-spin correlation function. The lines are fits to 
the power-law mentioned. 

defined by the Hamiltonian 

Hj = J^i-Sj, (27) 

where (ij) denotes a pair of nearest-neighbor sites. We 
will consider both the ID chain and the 2D square lat- 
tice (in both cases adopting periodic boundary conditions 
for systems an even number of spins). We also consider 
the J-Q model^ which includes four-spin interactions in 
addition to the exchange J: 

Hjq =Hj-QJ2 (S* ■ S, - ±)(S fc ■ S, - i). (28) 

(ijki) 

Here ij and kl form opposite edges on an elementary 
2x2 plaquette on the square lattice and the summa- 
tion includes both horizontal and vertical orientations of 
these edges on all plaquettes (i.e., the Hamiltonian obeys 
all the symmetries of the square lattice) . With the neg- 
ative prefactor of the Q term, this interaction clearly is 
related to enhancement of the parallel-dimer weight ci, 
in the notation of Fig. [TJ in the ground state wave func- 
tion (although the state is still of course not expected to 
be exactly reproduced by the CAP ansatz). 

A. Optimization method 

We start a variational calculation from some initial 
value of the parameters in the CAP state ([7]), typically 
a power-law form of the amplitudes h and all the cor- 
relation constants C = 1. When optimizing states for 



different values of some parameter (e.g., J/Q), we also 
normally start the calculation for a new parameter value 
from a previous calculation for some nearby value. One 
can also use this approach for different system sizes, al- 
though when increasing the system size initial values for 
the parameters corresponding to the longest bonds are of 
course not available and have to be set to some suitable 
values based on the longest previous bonds. In general, 
we have found that the starting point does not play an 
important role in optimization of AP and CAP states, in- 
dicating that the energy landscape is relatively smooth. 

To minimize the energy as a function of all parameters, 
we compute the energy and its derivatives. We then ap- 
ply either (a) the steeped-decent method or (b) a stochas- 
tic variant of it where only the signs of the derivatives are 
used, as discussed in Ref. [2(J A generic parameter p is 
in these two cases updated according to 

(a) p — > p — 5 ■ R ■ sign(dE / dp) , 

(b) p -> p- 5 ■ (dE/dp)/max(\dE\). (29) 

Here, in (a) R <G [0, 1) is a random number and in (b) 
max( | dE |) denotes the derivative that is the largest in 
magnitude among all the derivatives considered. The 
maximum shift 8 is gradually reduced, so that the vari- 
ational parameters eventually will converge. If <5 is re- 
duced sufficiently slowly, then one will reach a minimum. 
This minimum is not necessarily the global one, however. 
The stochastic scheme (a) should be better in avoiding 
local minimums, although one can of course never be 
completely guaranteed to find the global minimum. For 
the case at hand, the energy landscape appears to be 
relatively smooth, with no serious problems in consis- 
tently reaching the same minimum energy (regardless of 
the starting point, as mentioned above). Occasionally 
the method does fail, with independent runs of the same 
system leading to different final results. Typically, in a 
set of several runs, most of them will be consistent with 
each other, with only a small fraction of them deviating 
significantly from the majority value. As expected, the 
failure rate decreases when increasing the number of MC 
sweeps for sampling the VB configurations (leading to 
smaller error bars on the derivatives) and when reducing 
5 at a slower rate. 

We typically use a protocol based on an iteration num- 
ber k. For each k, the parameters are adjusted some 
number M of times (e.g., M = 100) based on derivatives 
obtained in MC simulations with oc k 2 steps. This way, 
the derivatives become more precisely determined as the 
solution approaches the minimum (where the derivatives 
decrease in magnitude, thus necessitating a larger num- 
ber of MC sampling steps to obtain statistically useful in- 
formation). The maximum parameter shift 6 in Eqs. (|29[) 
is of the form 5o/k a , with a = 3/4 a suitable exponent 
in practice, based on experience. 

Fig. [20] shows an example of the evolution of the energy 
and the sublattice magnetization squared in a run for a 
Heisenberg chain of length L — 256. It is clear that at the 
final step shown here, k = 30, the calculation has not yet 
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FIG. 21: (Color online) Staggered spin structure factor of 
the Heisenberg chain versus system size in variational AP 
and CAP calculations at different correlation levels. The re- 
sults are compared with exact results from unbiased QMC 
calculations.— 



FIG. 20: (Color online) Convergence of the energy per spin 
(top panel) and the staggered spin structure factor (bottom 
panel) as a function of the iteration number k for a ID Heisen- 
berg chain with L = 256 sites and all correlation amplitudes 
(i.e., r max = L/2 — 1) included in the CAP wave function (0. 
For each k, 100 adjustments of the parameters were carried 
out, each based on an MC simulation with 10 4 fc 2 updating 
sweeps. The points with error bars correspond to averages 
over the last 20 iterations for each k. The horizontal lines are 
results of unbiased QMC calculations.— 

completely converged (but note that one should not ex- 
pect convergence to the shown exact values, as the CAP 
states still are of course not sufficiently flexible to repro- 
duce the true ground state wave function completely). 
In practice, one can quite easily converge calculations for 
small systems essentially completely while large systems 
require very long runs. There is still some room for im- 
provement of the energy minimization protocol, as what 
we have explored so far are essentially schemes based 
on trial-and-error approaches. Nevertheless, the results 
to be presented next can be considered as almost op- 
timized and we do not anticipate that our conclusions 
would change based on more complete optimizations. 



B. Heisenberg chain 

An interesting question in one dimension is whether 
AP or CAP states can describe the critical ground state 
of the Heisenberg chain. As we saw in Sec. IIII1 with the 
simple parametrized ID AP wave function (|25p the cor- 
rect critical decay exponents (a = /3 = 1)$^ correspond- 
ing to this system cannot be achieved. In a variationally 
optimized state, we are not tied to any particular form 
of the amplitudes, and a sufficiently flexible variational 
wave function should then be able to capture the cor- 



rect criticality, including the logarithmic corrections that 
arise in the field-theory language due to a marginally ir- 
relevant operator. The question then is whether the AP 
or CAP states have this kind of flexibility, to possibly 
capture even such a subtle effect as the logarithmic cor- 
rections to the correlation functions.— 

To answer the above question, we have carried out en- 
ergy minimizations with the simple AP state (with all 
amplitudes as adjustable parameters) as well as with 
two types of CAP states. To include only the mini- 
mum amount of bond correlations beyond the AP state 
we include in ([7|) only the two ID bond-pair configura- 
tions with length-1 bonds; C{r\,r2) with r\ = ±1 and 
r2 = ±1. The case r± = l,r% = —1 corresponds to 
a single bond on a nearest-neighbor link and we can 
regard this as a normalization for the correlation fac- 
tors, (7(1, — 1) = 1, in the same way as we also set 
h(r = 1) = 1. There is then only one other corre- 
lation factor, C(— 1, 1) to optimize at this level. We 
also consider the extreme case of optimizing all C(ri, r%), 
r ma ,x = L/2 - 1, again with (7(1, -1) = 1. 

Let us first discuss the energy. As an example, for L = 
256 the exact energy per site is E/L = —0.44316 while 
for the AP state we obtained -0.44184. For the CAP the 
best energy when r max=1 is E/L — —0.44272, and with 
r max = L/2 — 1 it decreases to —0.44306. As discussed in 
the previous section, it is difficult to completely optimize 
long chains, so the optimal variational energies may still 
be somewhat lower. Following the trends as a function of 
system size, the relative energy error with the CAP state 
seems to remain at about 0.05%. 

Turning now to the spin correlations, in Fig. [5T] the 
staggered spin structure factor, defined according to 

L-l 

S-(tt) = £(-l) r C(r), (30) 

r=0 
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FIG. 23: (Color online) Levels of bond correlations. At level 
FIG. 22: (Color online) Optimized AP amplitudes for the n, the longest bonds (ri,r 2 ) for which the correlation weight 
Heisenberg chain of length L = 256. The line has slope -1.44. C(ri, r 2 ) in Eq. (0 is optimized (i.e., can be different from 1) 

are those marked by n. 



is graphed versus the system size for all the cases dis- 
cussed above. Since the exact C(r) decays as 1/r with 
a multiplicative logarithmic correction, the exact iS^tt) 
grows slightly faster with L than ln(L), as demonstrated 
with unbiased QMC results in Fig. [21] All the AP 
and CAP results exhibit a faster growth with L. When 
graphed on a log-log scale (instead of the lin-log scale 
used in Fig. l2"Tj) , the AP behavior is consistent with a 
power law, S(tt) ~ L 1_ " with a ps 0.70. With the CAP 
states, the data move closer to the exact points, but even 
with the maximally correlated CAP state the divergence 
is still somewhat too fast. 

It is also interesting to examine the optimized ampli- 
tudes of the AP state. Fig. [22] shows results for L = 256. 
Interestingly, a power law applies here for short and mod- 
erate bond lengths, with the deviations (enhancements) 
at large lengths likely related to the periodic boundary 
conditions (and some jaggedncss of the large-r data due 
to imperfect optimization, reflecting the total energy not 
being very sensitive to these "noise" features) . Even the 
r = 1 amplitude falls on the common power-law line in 
Fig. [22] i- e -: m the notation of Sec. IIIII the optimized 
state has A = 1. Looking at Fig. [6] when A = 1 the 
exponent a ps 0.75, quite close to a ps 0.70 obtained 
above with the optimized amplitudes. Thus, the bound- 
ary effects on h(r) seen in Fig. [22] appear to have only 
minor effects on the critical behavior. The conclusion 
for the optimized AP state is, thus, that a critical be- 
havior is reproduced, but with the wrong exponents for 
the correlation functions. Note, however, that a ps /3 for 
the applicable power-law obtained here, which is also the 
case for the true Heisenberg correlations (but with larger 
values, a = (3 = 1). 



C. Two dimensions 

We next systematically investigate the improvement of 
the energy with the inclusion of bond correlations in two 



dimensions, using several choices for the maximum bond- 
length r max in the correlation factors C(ri,r 2 ). Fig. [251 
illustrates all the bond shapes (r 1; r 2 ) at three correlation 
levels, with r max = 1, s/E, and 3 for correlation levels 1,2, 
and 3, respectively. 



1. Heisenberg model 

For the 2D Heisenberg model, previous variational AP 
calculations have shown that the energy error within this 
class of state is < 0.1% for large systems, and the spin 
correlations are reproduces to within 1% or better . 19 ! 20 
Although the system is strongly Neel-ordered and only 
has rapidly decaying short-range VBS correlations, in- 
cluding bond correlations with CAP states can still signif- 
icantly improve the energy further. Fig. 1241 shows results 
for L x L systems with L = 16, 32, and 64 at different 
correlation levels. The deviation from unbiased QMC 
calculations decreases with increasing correlation level. 
For L = 16 with r max = 3 the relative error is as small as 
ps 4 x 10~ 5 , while for the larger systems it is somewhat 
larger, about 10~ 4 . 

Going further and optimizing all correlation weights 
C(ri,r 2 ) with r < L/2 — 1, one should in principle be 
able to further improve the energy and obtain the best 
possible CAP state (with the kind of correlations in- 
cluded here) when L — > oo. The energy only improves 
marginally on the r max = 3 results, however. Fig. [25] 
shows results versus the system size for the energy as 
well as the sublattice magnetization. On the scale of 
the graphs, one can barely see any differences between 
the CAP and unbiased QMC results for L < 20, while 
for the larger systems there are some visible deviations. 
Here it should again be noted that the results for large 
systems are likely not completely optimized. As dis- 
cussed above in Sec. IV Al the energy depends only very 
weakly on the long-bond statistics, which implies that 
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FIG. 24: (Color online) Energy of the 2D Heisenberg model 
with variational CAP states at three different levels of bond 
correlations, according to the definition of the levels in Fig. 1231 
Level corresponds to the pure AP state, with no bond corre- 
lations included. The horizontal lines show energies obtained 
with unbiased QMC calculations (with the width of the lines 
corresponding approximately to the statistical errors). 
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MC evaluations of the corresponding derivatives are af- 
fected by relatively large fluctuations, leading to slow 
convergence. The sublattice magnetization is more sen- 
sitive to the long bonds, however, and this makes it very 
difficult to obtain completely unbiased results for large 
systems. For example, five independent optimizations 
for L = 32 with r max = 3 all gave the same energy 
within statistical errors, but the sublattice magnetization 
showed significant fluctuations, with the following re- 
sults: (m 2 ) = 0.1131, 0.1113, 0.1132, 0.1129, 0.1094, with 
the error bar approximately equal to 2 in the last digit. 
Here it can be noted that three of the results agree well, 
while two of them are clearly off. One may then conclude 
that the best optimized results should be around 0.1131, 
although to confirm this one should carry out a much 
larger number of independent runs. The correct results 
based on unbiased QMC calculations^ is (m 2 ) = 0.1128, 
less that 0.3% below the average of the above 3 consistent 
points. 



2. J-Q model 



FIG. 25: (Color online) Energy (top panel) and squared sub- 
lattice magnetization (bottom panel) of the 2D Heisenberg 
model obtained by unbiased QMC calculations and by opti- 
mized CAP state with all bond correlations included.— 
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FIG. 26: (Color online) Squared antiferromagnetic and VBS 
order parameters versus the correlation level in CAP states 
optimized for the Q model on a 32 x 32 lattice. 



As discussed in Sec. HI the J-Q model (|28|) exhibits 
a Neel-VBS transition at a critical value of the cou- 
pling ratio J/Q, with most precise estimate so far being 
(J/Q) c = 0.0447(2).^ An interesting question is whether 
this transition can be described by the CAP states. Here 
we consider the case J = 0, where the ground state 
is a columnar VBS. This VBS is very complex, how- 
ever, since the order parameter is only about 20% of 
the maximum possible value (i.e., for states with lcngth- 
1 singlets forming columns and no fluctuations around 



this configuration)^ The fluctuations are significant and 
have U(l) character up to a very large length scale (larger 
than what can currently be studied). As it turns out, the 
fully optimized CAP state (using r max = L/2 — 1) for this 
J = system does not reproduce the VBS order. Instead, 
as we will see below, the system is still on the Neel side 
of the quantum phase transition. 

First, let us again investigate the impact of including 
longer bonds in the correlation factor in Eq. ([7]). The en- 
ergy is improved very dramatically with increasing r max . 
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For example, for L = 32 the best optimized AP state 
has an energy E /N = -0.8013. The CAP states have 
E/N = -0.8215, -0.8229, and -0.8232, at correlation 
levels 1,2,3 in the scheme of Fig. [531 Going to r max = 5 
there is only a marginal energy improvement to —0.8233, 
which differs by about 0.1% from an unbiased QMC re- 
sult; E/N = -0.8240. 

Although an energy deviation of 0.1% would normally 
be considered excellent in variational calculations, the 
order parameters are still not well reproduced. Fig. [26] 
shows the dependence of both the Neel and VBS order 
parameters on r max . With increasing r max , the sublatticc 
magnetization is reduced and the VBS order parameter 
increases, as would be expected with CAP states in a 
VBS state. However, the VBS order parameter is still 
more than 30% too small at r max = 5, and it appears to 
be essentially converged at that point. Accordingly, the 
Neel order parameter is instead too large. 

Fig. [57] shows both order parameters calculated with 
r max = 3 as a function of the inverse system size, along 
with unbiased QMC results. Here one can observe that 
the agreement between the two calculations is very good 
for small systems, but the agreement becomes worse for 
increasing L. Asymptotically, the variational calcula- 
tions tend toward a weakly Neel ordered state, with no 
VBS long-range order, as opposed to the actual VBS 
ground state. This calculation serves to illustrate the 
insensitivity of the energy to long-distance correlations 
and the related difficulty in using the energy as a reliable 
measure of the quality of a state obtained by variational 
means (see Ref . [58| for a different example of this issue) . 

While this result could be seen as a failure of the vari- 
ational CAP states, it should be noted again that the 
actual Neel- VBS transition takes place at a very small 
coupling ratio, (J/Q) c ~ 0.045, and one cannot expect a 
variational calculation to reproduce the critical point ex- 
actly. For the CAPs considered here we have confirmed 
that the transition is pushed to a small negative J /Q, but 
we leave more detailed studies of the transition (which 
requires very well optimized states) for a future study. 



VI. SUMMARY AND DISCUSSION 

In conclusion, we have discussed VBS states and as- 
sociated quantum phase transitions in ID and 2D wave 
functions in the VB basis. VBS states appear naturally 
within the standard ID AP states, and we have here char- 
acterized the continuous Neel- VBS transition in such a 
class of states with amplitudes decaying as a power-law 
of the bond length. To stabilize a 2D VBS requires ex- 
plicit bond correlations beyond the AP states. We have 
introduced the CAP states, where correlations are en- 
forced through factors corresponding to bond pairs, as in 
the wave function ([Jj . We have shown how tuning of pa- 
rameters in 2D CAP states can lead to transitions from 
the standard Neel antiferromagnet to a VBS, in some 
cases with an intervening spin liquid phase. With the 




FIG. 27: (Color online) Squared antiferromagnetic and VBS 
order parameters versus the inverse system length for the 2D 
Q model within the CAP states with r- max = 3. 



parametrization considered here, the direct Neel-VBS 
transition is first-order, while the NeeHiquid and liquid- 
VBS transitions are continuous. 

The 2D Neel-liquid transition is of the same kind 
as in the pure AP states, although the short-distance 
VBS fluctuations are enhanced. Interestingly, the liquid- 
VBS transition is not associated with an emergent U(l) 
symmetry, although the columnar VBS (which is the 
VBS variant stabilized with the CAPs studied here) in 
principle supports this phenomenon^ and has been ob- 
served in QMC studies of J-Q models at the Neel-VBS 
transition . 33 ' 35 ' 52 Thus, in the states we have studied here 
the VBS should be induced by a relevant operator, in- 
stead of the dangerously irrelevant operator associated 
with the emergent U(l) symmetry in the "deconfmcd" 
criticality scenario. 

It remains an interesting challenge to find a 
parametrization of the CAP states such that a DQC 
point is obtained. In the study with short-bond correla- 
tions in Sec. lIVI we kept the single-bond amplitudes fixed 
with the form h(r) = 1/r 3 and varied only a correlation 
parameter p. In principle we could also use h(r) = l/r K 
and also vary k. It would be interesting to study the full 
phase diagram in the plane (p, n) for the two parametriza- 
tions of the bond correlations, and also with different 
choises of the correlation factors. 

We note a recent study of AP states, with a cer- 
tain parametrization of the amplitudes, to describe the 
Neel to quantum-paramagnetic transition in the bilayer 
Heisenberg model^ The correct 3D classical Hciscnberg 
exponents were obtained for this transition, i.e., the AP 
states contain effectively long-range interactions that al- 
low (2+l)-dimensional criticality to be reproduced within 
a 2D configuration space. In principle it seems that this 
should be possible to achieve also for the DQC transition, 
and if such a program to construct a simple CAP wave 
function is successful, it would likely lead to further use- 
ful insights into the mechanism of spinon deconfinement 
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and emergence of the effective U(l) gauge field. 

It is possible that CAP states with all parameters ad- 
justed to minimize the energy of a model such as the 
J-Q model could lead to a DQC point. Here we car- 
ried out such variational QMC calculations for the stan- 
dard J-Q model with four-spin interactions. This model 
has a critical point very close to J/Q — 0, however, 
and in the variational calculation the J = system is 
still inside the Neel phase. In principle one can still 
study the phase transition by going to negative J, but 
in that case Marshall's sign rule cannot be proven rig- 
orously (although most likely it should still hold when 
J/Q is small and negative). We will investigate this case 
in future studies, and also consider the model with six- 
spin interactions,— where the transition point is at much 
larger J/Q and should remain well within the range of 
positive J/Q values within optimized CAP states. If in- 
deed the correct type of criticality can be achieved, then 
by examining the bond amplitudes and correlation fac- 
tors in the optimized state it may also be possible to 
construct a class of CAP states with a single tunable 
parameter (instead of all the parameters changing as a 
function of J / Q in the variationally optimized states) to 
drive this type of criticality — which the paramctrization 
used in the present paper was not capable of. 

Beyond the case (I) and (II) parametrizations of the 
CAP states that we considered here, we have also ex- 
plored other cases. In particular, with c\ = C3 = p and 
C2 = C4 = 1 /p, in the notation of Fig. [TJ we have found a 
plaquette VBS, in contrast to the columnar VBS states 
obtaining in cases (I) and (II). In future studies it will 
also be interesting to study the nature of the transition 
from a Neel antifcrromagnct to a VBS in this case. 

It is in principle possible to further improve on the 
CAP states, by including more correlations. We here con- 
sidered pairs of bonds connected to a nearest-neighbor 
link only. We have verified that this gives a better 
variational energy (for the Heisenberg and J-Q mod- 
els) than next-nearest-neighbor links, and therefore, most 
likely, the nearest-neighbor links are optimal for intro- 
ducing correlations in this kind of CAP states. One could 
also combine several types of correlation factors, and in- 



clude also factors for correlations between more than two 
bonds. In practice this may not be worth the effort, how- 
ever, as the main utility of CAP states should be (i) to 
have a simple class of states to capture the Neel- VBS 
transition and (ii) to use them as "trial states" for pro- 
jector QMC calculations in the VB basis .-i^^ While 
the variational states can be improved, in practice it is 
better to project out the ground state exactly using QMC 
if completely unbiased results are needed, and too many 
parameters in a CAP state defeats the purpose of (i). 

A very interesting question is whether the phase tran- 
sitions we have discussed here can be realized in ground 
states of reasonable Hamiltonians, with only local in- 
teractions. In particular, the Neel-liquid-VBS series of 
phases would be of interest in this regard. We already 
know from the worh of Cano and Fendley on the short- 
bond RVB that there is a local parent Hamiltonian for 
that state. It is then also appears very plausible that 
some local pertutbations of this Hamiltonian will effect 
the stiffness constant characterizing the RV B 21 i 55 and 
governing its critical VBS correlations. Thus, a class of 
local Hamiltonians should be able to capture the whole 
spin liquid phase in our case (II). Then, it also seems 
plausible that other local perturbations can drive the sys- 
tem into a Neel or VBS states, e.g., the Q term of the J-Q 
model should do this. Since the Cano-Fendley Hamilto- 
nian has a sign problem in QMC calculations, some other 
methods would be needed to study phase transitions in 
perturbations of it. 
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